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The Doi formalism treats a reaction-diffusion process as a quantum many-body problem. We use 
this second quantized formulation as a starting point to derive a numerical scheme for simulating 
A' — > reaction-diffusion processes, following a well-established time discretization procedure. In 
the case of a reaction zone localized in the configuration space, this formulation provides also a 
systematic way of designing an optimized, multiple time step algorithm, spending most of the 
computation time to sample the configurations where the reaction is likely to occur. 

1. INTRODUCTION 

The reaction diffusion process X — > with space dependent reaction rate, models a number of interesting chemical 
and biological situations. It models for instance the capture dynamics of a mobile ligand diffusing in a configuration 
space until being captured by a static receptor. This configuration space may not just be the three dimensional 
volume containing the molecules which participate to the reaction, but must be considered in a wider perspective. 
For example, in the case of a ligand bound to one end of a polymer spacer while the other end is grafted onto 
a surface (bridging situation [![), the diffusive dynamics naturally takes place in a configuration space taking into 
account all the internal degrees of freedom of the chain, even though the reaction rate mainly depends on the position 
of the polymer end bearing the ligand 0- Biology and bio-technologies provide numerous examples of adhesion 
mechanisms with kinetics limited by similar reaction-diffusion dynamics [H, Q . The most standard immunodiagnosis 
assays utilize the aggregation processes based on antigen-antibody specific adhesion. It is now possible to set up finely 
tuned experiments with a high enough temporal resolution for the comparison of reaction kinetics with analytical or 
numerical predictions @ . 

Polymer chemistry also offers a wide set of situations which can be advantageously modeled by a X — > process. 
The cyclization reactions occurring during macromolecules synthesis are the examples 0, H| ■ The situation can be 
further complicated if an external hydrodynamic flow is present, or more generally, if conservative or dissi pat ive force 
fields are known to act upon the system under study , then calling for intensive numerical simulations jlOl . ITl| . 

We present in this paper a general method for finding efficient numerical simulation algorithms, using a non 
hermitic Hamiltonian formulation of the stochastic process introduced by Doi [13 ]. and based on a Trotter splitting 
of the discrete time evolution operator. This strategy has proved successful when applied to Brownian and DPD 
(dissipative particle dynamics) dynamics, leading to the generation of algorithms, with controlled time-reversal or 
accuracy properties |l3l4lq . 

Moreover, we suggest for this problem an original multiple time step algorithm, based on some internal state 
switching A •<-> B, which comes as a straightforward extension of the physical reaction-diffusion process. We then use 
a different integration time step, depending on the internal variable state, either fast or slow. 

The essential strategy is to introduce fast coarse grained sampling away from the reaction zone (while keeping slow 
sampling near the reaction zone) in order to achieve the desired computational efficiency and accuracy. This algorithm 
comes as a new member in a family of optimized multiple time step algorithms, such as the RESPA [l7l . Ila ]. and 
previously restricted to a fixed number of moving particles. 

In this work, we restrict our approach to the case of X — > non-interacting particles. However it can be easily 
extended to more general cases without difficulty. Although the formalism of Doi is not new, it is, to our knowledge, the 
first time that it serves as a starting point for designing a specific integration scheme for reaction-diffusion simulations. 

We introduce in Section 2 the Doi formalism [12|, |l9( and derive the corresponding single time step algorithm. We 
underline in Section 3 the danger of naive incorporation of multiple integration time-steps for different regions of the 
configuration space, and show how to properly formulate multiple time step algorithms. 

Section 4 demonstrates the validity of the whole scheme for particles diffusing on a spherical surface. The simulation 
time can be reduced by a large factor without noticeable artifact by choosing proper step sizes. The estimate of the 
expected speeding up efficiency is provided. Perspectives arc outlined in our conclusive Section. 
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2. FORMALISM 



The standard master equation for the X — > process is 

<9 t n(r, t) = £>V 2 n(r, t) - fiV(n(r, t)f(r, t)) - S(r)n(r, t), 



(1) 



where n(r,t) is the probability of finding a particle at a position r, at time t, S(r) is the space dependent reaction 
rate, D the diffusion constant, fj, the mobility coefficient, linked to D by the Einstein relation D = ksTfi, and f(r) a 
space dependent, constant force field. Eq. ([IJ is a Fokkcr-Planck-Smoluchovski diffusion process, with an additional 
sink term accounting for particle decay. Usually, the rate S(r) is a sharp function of r, and represents, for instance, 
the very localized region in configuration space where either a binding or a chemical reaction can take place, with a 
precision lying in the range of the angstrom. The diffusive behavior underlying cq. ([T]) corresponds to an overdamped 
Brownian diffusion, in which the inertia of the diffusive particle plays no part. 

The X — > process does not conserve the number of particles, but is not a many-body problem, as there is no 
interaction affecting either the spatial motion or the reaction rate. Thus, following the idea of Doi, we introduce the 
second-quantized time evolution operator, that we subsequently call "Hamiltonian" by reference to standard quantum 
mechanics [HI, Qj], [2(J (it is also referred as a Liouvillian in related ref . [2l[ ) . 



V.= dr 



D(Va\r) • Va(r)) - fx{f{r) ■ Va^r))a(r) + S'(r)(a t (r)a(r) - a(r)) 



(2) 



In eq. @, the first term corresponds to the Brownian diffusion term (Hd), the second one to the force induced drift 
term (Hf), and the last one is the space dependent reaction X — > (Hr). The operators a(r) and (r(r) are bosonic 
creation and annihilation operators. Here, the use of bosonic operators means that multiple occupancy at a given 
position r is possible. However, in the absence of reverse reaction — > X, or any other reaction creating particles X, 
the total number of particles present in our system can only decrease. If one starts at t = with a single particle at 
the origin r = r (i.e. n(r 7 r ) = S(r — r )). 



|$(0))= O T(r )|0) 



(3) 



then the evolution operator exp(— T-Lt) of our system, either moves or destroys this particle, and double occupancy 
never happens. Thus, for a single particle reaction dynamics, formally it makes no difference here whether one chooses 
to use bosons or fermions in eq. ((2|). 

The connection between the second-quantized formalism and the standard formalism is made with 



n(r,t) 



|exp 



dr'a(r'))a(r) |$(t)} , 



(4) 



where = exp(— Ht) |$(0)) represents the state of the system at time t [2(j. As double occupancy is precluded, 

one can safely replace the above expression by 



n(r, t) = (0| a(r) |$(t)) = (0| a(r) e -^ 1 |$(0)) 
and define a propagator, for positive values of t 

G(r 2 ,rx;t) = (0| a(r 2 ) e - ?i *a t (r 1 ) |0) . 



(5) 



(0) 



In eq. ([2]), the Hamiltonian is a sum of three non-commuting terms Jio, T~Lf and %r, and exp(— Tit) is not known 
analytically. Following the Trotter-Strang idea, we split the evolution operators as 



-{H D +n F +H R )t 



t/At 



t/At 



e 2 



Hi 



-AtH L 



e 2 "-f . e 2 



Hi 



(7) 



The approximation becomes exact for At — > 0, and the above splitting is known to be exact up to terms of order 
At 3 . Expression ([7]) has a clear algorithmic interpretation. It tells us to iterate a large number of times (N = t/At) a 
sequence of five steps, which altogether are supposed to reproduce faithfully the dynamics of the process over a time 
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interval At. The error committed when doing this approximation is of order At 3 and originates from the fact that 
we have split the exponential of a sum of three non-commuting terms. 

We now consider each one of these steps: exp(— A£Hr/2) represents the action of the reaction term, at a fixed 
position, during a time interval At/2. Then, the survival probability of a particle initially at r is exp(— S(r)At/2), 
and the probability of having reacted is just 1 — exp(— S(r)At/2). 

The term exp(— AtT-^/2) represents a pure drift along the (constant in time) force field f(r). This is tantamount to 
solving the first order differential equation r(t') = /i/(r(t')), with initial condition r(t' = 0) = r, during a time interval 
At/2. We can check that, at the desired order of accuracy (i.e. At 2 ), this corresponds to a shift r — >• r + fj,Atf(r)/2, 
where the pseudo force field f(r) is closely related to the original force field f(r), according to the Heun, or Collatz 
scheme (see Appendix II and ref. [13 )■ 

The term exp(— AtW-n) represents a pure Brownian diffusion in flat space, which is nothing but a Gaussian (Wiener) 
process with variance 2D At. 

The splitting in Eq. (JT]) can be cast in the algorithm presented in Table fl] At the end of this procedure, the position 
of the particle is 

+ Ar D , (8) 

and the survival probability is exp(— [S(r(t)) + S(r(t + At))]At/2). If a reaction takes place, the reaction time must 
be set to 

At 

tfi = t+— , (9) 

and the reaction diffusion process stops. In fact, we cannot define the reaction time more accurately than that, 
because any attempt to link the reaction time to one particular step in the previous algorithm would be pointless. 
The five steps of Algorithm I are the consequence of a formal, mathematical splitting of ([7]), and does not represent 
a consecutive sequence of events in real time. 

Such an algorithm should prove successful when applied to smooth, diffcrentiablc space varying functions S(r) and 
f(r). The accuracy of its predictions is expected, when averaged over the different Brownian paths and random 
reaction times tn, to reproduce the results of the exact, continuous time process, up to and including the order At 2 . 
These algorithms are known, in applied mathematics, as "weak order two algorithms" [23] ]. 

3. OPTIMIZATION AND MULTIPLE TIME STEP ALGORITHMS 

The accuracy of the above numerical scheme is limited by the regularity of the functions S(r) and / (r) appearing 
in the problem. The characteristic length scales As or Af upon which these functions vary set the upper bound for 
the discrete spatial increment Ar, which in turn set the upper bound for the discrete time step At of our dynamics. 

In capture problems, it commonly happens that the reaction zone occupies only a tiny fraction of the configuration 
space, with a length scale As of the size of the angstrom. 

Using a large time step At enables a fast sampling of the configuration space. However the risk of missing the 
reaction zone increases accordingly overlooking the possibility for the Brownian path, which is a continuous process, 
to come across the reaction zone and react (Fig. []} . 

In practice, the time step At must be chosen in order to keep the typical length step | \r (t + At) — r (t) \ | smaller than 
As by a factor two or three. On the other hand, the choice of a small time step can make the computational effort 
unbearable if L 3> As, and when nothing much happens outside the reaction zone. For those reasons, it is desirable 
to have a multiple time step algorithm at our disposal. 

A naive attempt to make a multiple time step algorithm is to divide the configuration space into two zones: a slow 
zone where one applies the algorithm with a small time step At, and a fast zone where one applies the algorithm with 
a larger time step J- At, where T is an arbitrary acceleration factor. Convenient choices for T are powers of 2, such 
as J- = 8, 16 . . ., to enforce commensurability between small and fast dynamics. 

As a test case, we apply this idea to a one dimensional particle diffusing in the interval x G [—0.5, 0.5] (periodic 
boundary conditions) with "slow" and "fast" domains and an acceleration factor T = 8. The large time step applies 
for a particle in \x\ > 0.2 and the small step for particles in \x\ < 0.2, which would correspond to a reaction zone close 
to x = 0. Figure [2] shows the result of this simulation. A sharp discontinuity between fast and slow zones leads to 
an inhomogencous density of particles. This artifact is located near the slow/fast boundary, with width ~ \Z2TDAt 
given by a diffusion length, and due to a local breakdown of the detailed balance relations, where incoming particles 
large moves are not balanced by the reverse small moves. 



r(t + At) =r(t) + 



fiAt 



f(r(t)) + f(r(t) + fif(r(t))At/2 + Ar u ) 
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In order to suppress this artifact, we propose to introduce a region, called "exchange zone" , of size A w where fast 
and slow dynamics overlap. This is achieved by adding to the diffusing particles an internal state variable, a = A or 
B, switching randomly between a state A and a state B, like in a reversible chemical reaction A <H> B. 

The standard master equation of the process becomes: 

d t n A (r,t) = DV 2 n A (r,t)-iJ,V(f(r)n A (r,t))-S(r)n A (r,t)) 

+w AB (r)n B (r, t) - w BA (r)n A (r, t)\ 
d t n B {r,t) = DV 2 n B (r,t) - fiV(f(r)n B (r,t)) - S(r)n B (r,t)) 

-w AB (r)n B (r,t) + w BA (r)n A (r 7 t). (10) 

The internal state (cr = A or a = B) does not alter the diffusion nor the reaction properties of the particle, a = A 
or B is a hidden variable, while all the space dependent properties remain unchanged compared with eq. JT]). We 
loosely call particle A a particle with a = A and particle B a particle with a = B. In the slow diffusion limit, the 
equilibrium molar fraction of particles A and B is fixed by the exchange rates w AB (r) and w BA (r), and given by 



x A 



x B 



»A 



w AB {r) 



n A + n B 
n B 



w AB (r) + w BA (r) 
w BA {r) 



n A + n B w AB (r) + w BA (r)' 



(11) 



The exchange rates w AB (r) and w BA (r) can be arbitrarily chosen, and tailored to compromise between an strong 
conversion rate (w must be large), and a relative smooth spatial variation (w cannot vary too sharply with r). We 
keep the size of the exchange zone A w two to three times larger than the spatial increment of the particles to preserve 
the total particle number. 

The second-quantized Hamiltonian corresponding to this situation is: 



% a = % A M, + 'K A ,D + W-A.F + 'H-A.Ex + 'HbM. + "Hfi.D + ~Hb,F + T~L Bt Exi 



with components 



T~La,r 
Ha,d 
Ua,f 
Ha, Ex 

7~Lb,d 
Hb,f 

Wb,Ex 



dr 
dr 
dr 
dr 
dr 
dr 
dr 
dr 



S(r)(a^(r)a(r) - a(r)) 
D(Va f (r) • Va(r)) 
-li(f(r)-Vat(r))a(r 
WBA{r) ^a^(r)a(r) — b\r)a(r) 
S(r)(tf(r)b(r)-b(r)) 
D(W(r)-Vb(r)) 
- M (/(r)-V6 t (r))6(r 
wab{i") I 6 t (r)6(r) — a)(r)b{r 



(12) 

(13) 
(14) 
(15) 
(16) 
(17) 
(18) 
(19) 
(20) 



where b, b^ are bosonic creation/annihilation operators associated with particle B, and a, associated with particle A. 
Up to this point, particles A and B obey the same dynamics, and one is not faster than the other. The numerical 
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scheme originates from the splitting of exp(— F1-L a Ai) as follows. 



exp(-JTCAt) 



x e 




{ 



} 



■772 



x 



e -At-H B , Ex /2 e -AtW B ,R/2 e -AfHB,jr/2 e -At«B,D e -AtH B ,jr/2 e -At«B,B/2 e -AtH B ,Ex/2 



(21) 



The above splitting has a straightforward algorithmic interpretation. These terms indexed by B only acts on the 
particles whose internal state is B and leave the A particles unchanged with both position and internal states. The 
converse is true. 

One recognizes three occurrences of the splitting (JT]), enclosed in brackets [•], which do not depend on the internal 
state of the particles, but only on the time interval (respectively At or J- At). The splitting procedure has already 
been detailed in Section [5J Therefore, any implementation of the single time step algorithm can be reused in the 
multiple time step situation, while providing also a natural checking procedure for this more complex algorithm. 

The effect of cxp(— AfHB.Ex) is to exchange the internal state of B particles to A, with a space dependent probability 
equal to 1 — exp(— Aiuus (?")), corresponding to a waiting time At. Meanwhile, particles A are not affected and stay 
at the same position r. The resulting algorithm is summed up in Table [TT1 

The dynamics of particles B is treated with a discrete time step At, while the dynamics of particles A is treated 
with a larger time step TAt. It is tempting to call A the fast particle (or the fast state) and B the slow particle (or 
the slow state). This linguistic shortcut must be used with care, because A and B particles share identical dynamic 
properties, and this is only the way we treat them in our discretization procedure which differs. 

In practice, the exchange rates must be tailored in such a way that the reaction zone falls into a region where 
particles B dominate. If we arc confident that particles with state A never wander around the reaction zone, we can 
safely drop the terms exp(— / HA,nAt/2) from the algorithm. Moreover, there is a considerable room for optimizing 
the multiple time step algorithm, at the expense of a slightly more complex programming. 

If necessary, the idea can be further extended by adding more internal states A\,Ai, A3..., corresponding to 
consecutive nested sub-domains of the configuration space. The corresponding generalization of the algorithm is 
straightforward. 

To check the relevance of the above idea, we now introduce an exchange zone for particles diffusing in one dimensional 
space as described earlier. The "slow" particles and "fast" particles exchange their status in a slab 0.2 < |x| < 0.3 
of width A w =0.1. The acceleration factor is kept the same T = 8. The average density of particles A and B now 
overlaps (see Fig. [3]), and the artifact shown in Fig. [5] seems to be removed. Thus, our splitting (|2ip looks like an 
efficient strategy for keeping the advantages of a multiple time step algorithm, without suffering from appreciable 
unphysical features. 

The multiple time step algorithm shows equal efficiency in higher dimensions. To demonstrate the advantages of 
the exchange zone in higher configurational space, we consider independent particles diffusing in a d s dimensional 
configuration space. Particles are reflected by a (hyper) spherical outer boundary located at distance r = R e from 
the origin and become captured when hitting an absorbing (hyper)spherical inner boundary at r ~ Ri, as depicted 
in Fig. 3J This system is a straightforward generalization of the previous unidimensional case. Despite inherent 
complexity, this reaction-diffusion situation remains amenable to an analytical exact solution thanks to the rotational 
invariance of the system. Details of the analytical solution are provided in Appendix [Cj 

For the multiple time step procedure, we use T — 8 and set two exchange zones of width A„, = 0.15 at the vicinity 
of the inner and outer walls (respectively Ri — 1. and R e ~ 2.5). Figure [5] shows no visible artifact in the particle 
distribution when the exchange zone width A^, is kept finite, while the distribution is no longer uniform if A w = 
as shown in Fig. [6l In Fig. we compare the numerical and analytical Laplace transforms of the survival rate 
PD,Ri,Re(s) (c.f. Appendix . starting from a uniform particle distribution at initial time, for d s = 3, 5, 8 and 10. 
The combination of an exchange zone and internal switching state dynamics maintains a quasi-uniform distribution, 
while reducing the computation time by factor 2. This fits our expectation as discussed in the next section. 



4. 



EXAMPLE 



We now turn to a realistic situation to see whether our strategy keeps its promises. We consider, for example, two 
localized colloid carrying a ligand and many receptors, respectively, on their surfaces. These colloids combine in the 
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form of doublet upon ligand-rcceptor binding. Wc model the rotational diffusion of colloids as random walks of a 
particle on a sphere. The motion of a Brownian particle evolving on a sphere is parameterized in polar coordinates 
(9, <fi), with a capture zone defined by 9 < 9 C (Fig. HJ. We performed a statistics of survival times t s , starting from a 
uniform distribution of initial conditions over the sphere, until the particle encounters the capture circle 9 = 9 C . An 
analytical solution of this problem, presented in Q served as a benchmark for checking out the whole procedure. 

To implement the optimization procedure, we defined an exchange zone between angles 6i n j = 0.4 and 9 sup = 0.6, 
while the capture zone was set to 9 C = 0.03. Figure |H] shows the survival probability as a function of time for various 
acceleration factors J- = 4, 8, 16 and 32 without any appreciable deviation, while the angular distribution of A and B 
particles is presented in Fig. [TO] Meanwhile, we must check whether the multiple time step algorithm preserves the 
uniform distribution of positions when there is no capture (reaction off), which was achieved by counting the number 
of recurrences, as shown in Fig. 1111 

For a proper choice of we must keep the length step smaller than the width of the exchange zone A w . Choosing a 
larger T might accelerate the fast particle but requires an accordingly larger exchange zone, where the particle might 
spend a longer time. Hence, we observe that the CPU time is reduced by factor 6 for T = 16 and by 9 for T = 32 
with A w = 0.2 and 0.4 respectively. Corresponding CPU times are provided in Table IIIII For a given J 7 , the time 
step of fast particle is increased by J 7 , and the corresponding length step (A9) by \[T . 

A rough estimate of the CPU time dependence on T can be obtained with a few approximations. Provided that 
the volume of the interior zone (including exchange zone) where particles are "slow" occupies a fraction a of the total 
configuration space volume L d , and that the remaining part of the configuration space is filled with "fast" particles 
subject to the acceleration factor J 7 , then the time needed to sample the whole configuration space is approximately 
proportional to 

1 — a , 
a+—, (22) 

and the computing time reduced by a factor J 7 / (1 — a+ To). In our case, with a the sphere of area 4-7T and an exchange 
zone [0.4,0.6], we observe a confirmation of the expected behavior (see Fig. IT2"j) . Equation (|22|) also accounts fairly 
for the CPU time associated with the diffusion in d s dimension illustrated in Fig. U with a coefficient a comprised 
between 0.3 and 0.6. 

If one now wishes to compare the CPU time of the multiple time step algorithm with exchange procedure (Table |T| 
with its single time step counterpart (Table [TTJ) , eq. (f2"2"j) must be changed in 

(1+e), (23) 



T 



where e accounts for the relative increase of computations associated with the exchange procedure, which depends on 
each specific case but should remain small in all relevant situations. 

The acceleration factor T cannot exceed an upper value T max dictated by the exchange zone size A„ and slow zone 
size La 1 ' 11 . Thus, increasing significantly J- without reintroducing spurious sampling features requires the enlargement 
of both exchange and slow zones, and hence does not improve so much the efficiency. The optimal strategy may consist 
in partitioning the configuration space into nested domains A\ , A2 , A3, . . . associated respectively to acceleration factors 
1, J- 2 . . . and with geometrically increasing characteristic sizes. A well balanced partition should be such that the 
computer time linked to crossing each domain Ai remains approximately constant (a procedure reminiscent from the 
so-called Monte-Carlo umbrella sampling [l4[ or Transition Interface Sampling methods |24|). In practice, values of 
T equal to 8 or 16 sound promising. 

Finally, we considered the case of Nl ligands simultaneously diffusing on a sphere covered with Nr randomly 
placed receptors. There are Nr capture zones and each capture zone is defined as angular distances from each 
receptor. As discussed in ref. [6j, one expects the survival probability of the ligands to decay exponentially at long 
times, with an inverse relaxation time proportional to the number of receptors. In Fig. 1131 wc compare the survival 
time distributions obtained with the single time step and multiple time step procedure. As shown in Fig. 1131 this 
distribution is insensitive to the choice of the simulation procedure with the choice of a proper J- (here J- = 16). 

The computational time of three different procedures is compared: single time step algorithm, multiple time step 
algorithm with exchange zone and multiple time step algorithm without exchange zone (naive procedure). The 
benefits of the multiple time step algorithm decrease when the number of ligands Nl increases. When a single 
receptor (Nr = 1) is present on the sphere, the multiple time step algorithm performs always better than the single 
time step one (Fig. ITU . If more receptors are present (Nr = 10), the benefits of the multiple time step procedure 
with exchange zone are lost when the number of ligands reaches Nl — 20 (Fig. [T5| . This is because the switching of 
the particle status in the exchange zone requires as many as Nl x Nr operations, which cancels out the gain from 
using the acceleration factor T x Nl (large e in expression (|23[) ). However, if pairwise interactions among ligands 
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(receptors) were to be considered, the additional number of operations for particle status switching would become 
negligible compared with the quadratic iV£ number of operations required by the two-body interactions (small e in 
expression (|23|l ). Thus, multiple time step algorithm is again clearly favorable. 

5. CONCLUSIONS 

We presented a systematic procedure to simulate the Brownian diffusion of independent (non-interacting) and non 
conserved particles. Based on the second-quantized formulation of Doi, our approach uses a Trotter splitting technique 
and provides a weak order two algorithm for simulating the original process. 

We also showed how this reaction-diffusion approach could help in finding an optimized and faithful sampling of 
the configuration space, in cases where the spatial resolution of the Brownian paths is not necessarily homogeneous, 
without showing the artifacts that other naive approaches encounter. 

We tested our model on particles diffusing on a sphere and reacting at the vicinity of a pole. We were able to 
cut by an order of magnitude the computing time, while preserving the accuracy of the numerical scheme. When 
many particles diffuse simultaneously in the presence of multiple reacting zones, the efficiency of the multiple time 
step procedure declines for large number of particles. Although the efficiency depends on the overall algorithmic 
complexity of the various parts of the simulation scheme, the computation time of multiple time step algorithm is 
bounded by the computation time of single time step algorithm. 

In the multiple time step algorithm, our trick consists in introducing an extra state variable which stochastically 
switches between two internal states. This procedure is certainly not restricted to the Brownian dynamics exemplified 
above, but might be also of interest for deterministic, or inertial Langevin dynamics. This formalism will be the 
natural starting point for such generalization. 

Acknowledgement Author F.T. thanks KIAS for financial support during its visit. N.-K. L. acknowledges financial 
support by the Korea Research Foundation Grant funded by the Korean Government (MOEHRD, Basic Research 
Promotion Fund) (KRF-2005-204-C000024). 
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Appendix A: Connection between second-quantized formulation and master equation 

Starting from 

n(r,t) = (0|o(r)e-**|$(0)} (Al) 

one splits the evolution operator 

e -nt = e -H(t-t'). e -jit' (A2) 
and introduce a projector (valid only in the subspace of states with no multiple occupancy) 

1 = ( |0) (0| + J cW(r) |0) (0| a(r)^j , (A3) 

leading to 

n(r, t)= J dr' (0| a(r)e-^*-*V(r') |0) (0| o(r / )e - **' |*(0)) (A4) 

and to the well known propagator structure of the time evolution : 

n(r,t) = J dr'G(r,r';t-t')n(r',t') (A5) 

In the same way, the partial differential equation (TTJ) is obtained by differentiating eq. (|A1[) : 

d t n{r,t) = - (0\a(r)'He- ilt |$(0)) 

= -(0\Ha(r)e-™\^0)) + (0\[H,a(r)]e-™\$(0)). (A6) 

The Hamiltonian is a sum of normal ordered (creator operators on the left of annihilator operators) terms and 
pure annihilator terms a(r), which ensures that only the second term of (|A6|) is non zero. A standard commutator 
calculation leads to cq. (fTJ). 

The connection between the second-quantized formulation of the dynamics and the algorithm is done by splitting 
the evolution operator, such as in (J7J or in cq. (|2I[) . and by introducing as many projection operators (|A3I) as 

necessary between the exponentials. Each one of the (0| a{r{)e~ HAt a) '(1*2) |0) has a straightforward interpretation, 
and corresponds to an elementary step of our algorithms I and II. 



Appendix B: Heun and Collatz schemes for a first order differential equation 

For integrating an ordinary differential equation (o.d.e) r = f(r), the simplest numerical scheme is the Euler 
method : 



r{t + At) = r{t) + Atf(r(t)). 



(Bl) 



However, this scheme is known to be an "order one" integrator, which in the present situation, means that the 
approximated solution departs from the exact one, on a time interval At, by an amount of order CAt 2 . Clearly, this 
error is not acceptable as it would spoil the accuracy of the whole algorithms I and II. 
Among the number of existing "order two" integrator, we propose to use either : 



, . At 
r(t + At) = r(t) + — 



f(r(t)) + f(r(t) + Atf(r(t)) 



(B2) 



(Heun scheme), or 



r{t + At) =r{t) + At 



(B3) 
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(Collatz scheme) [22|]. These schemes both require two evaluations of the function f(r) each, per time step. Once 
At has been chosen, it is possible to give to eq. (|B2[) and (|B3|) a form similar to eq. (|B1|) . by defining a pseudo force 
field : 



/» = \ 



/(*■) + / r + /(r)At 



(B4) 



(Heun pseudo force field), and : 



(B5) 



(Collatz pseudo force field), such that 



r(t + At) = r{t)+Atf(r(t)). 



(B6) 



becomes a genuine order two integrator. We suggest to use either (|B5[) or (|B6I) in the steps of the algorithms associated 
to e-^ A *. 



Appendix C: Diffusion and absorption in d s dimensions 



To compute the survival probability between an absorbing boundary at r = Ri and a reflecting boundary at r = R e , 
we introduce the concentration of particle c(r, t) at position r and time i, obeying: 



(d t - A)c = 0, 



with c(Ri,t) = 0, d r c(R e ,t) = and where 



A 



1 d 



r dr dr 2 



(CI) 



(C2) 



is the hyperspherical Laplacian in d s dimensions, discarding everything but the radial variable r. In this example, the 
diffusion constant D is set to one without loss of generality. Denoting cq = c(r, 0) the initial uniform concentration, 
the time Laplace transformed c(r, s) obeys 



(• 



1 d 



r dr dr 2 
2w d 3 /2 -R 



)c(r,s) = c ; 
co = 1, 



drr rf 3 -i 



r(d./2) Jr, 

the solution of which can be expressed in terms of modified Bcssel functions, with v = d s /2 — 1: 



c(r, s) = — 
s 



r" /„+i(Vs J R e )^(Vs-Ri) + ^+i(^ e )^(A/i-Ri) 



(C3) 
(C4) 

(C5) 



Integrating c(r, s) over r between Ri and R e gives the Laplace transform Pd a ,Ri,R s (s) of the survival rate Pd 3 ,Ri,R e (t): 



Pd B ,Ri,R e (s) 



IiU {yfsRi)K ^{sfsR e ) - {yfsR e )K^ {^fsR l 



\fsRi 



I ia .(y/8R e )K ix _ i {y/8R i ) + K* a .(y/8R e )I* a ._ 1 (y/SR i ) 



(C6) 



The survival rate Pd a ,Ri,R e (t) is a monotonically decreasing function of time t, giving the proportion of particles 
remaining after t, starting from a uniformly distributed random initial position in configuration space. It is related 
to the average survival time (t) by the simple relations: 



(*) 



lo 



&tPd.,R i ,R e (t) = lim Pd,,R,,Re(s). 
s— f0 



(C7) 
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The simulation provides us with a finite sample of N survival times {U}, 1 < i < N. These survival times can be 
reordered into a monotonically increasing sequence ta-u) by means of a suitable permutation a of the indexes. Such a 
numerical sample approximates the survival rate Pd s ,Ri,R s (t) up to statistical fluctuations of order 1/y/N. 

In practice, this is achieved by plotting 1 — i/N as a function of Alternatively this sample can be used to 

numerically obtain the Laplace transform of the survival rate Pd s ,Ri,R c (s), by plotting p( num )(s) versus s, as shown 
in Fig.0 with p( num )( s ) defined as 

p(nu m)(s)= l£l_CXp(- i ^ (cg) 



Figures and Tables 
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Algorithm I: reaction-diffusion over At 

1. Consider the probability that a particle located at r(t) reacts during a time interval At/2, 
given that the probability of this event is 1 — exp( — S(r(t))At/2). If the particle reacts, end the 
procedure and set the reaction time tR to t + At/2. 

2. Shift the position of the particle by /i/(r(f)) At/2, as a consequence of the drift force. 

3. Draw from a Gaussian distribution of variance 2D At a random step Ar D , and shift the 
position by Ar D . 

4. Shift the position by nf(r(t) + fif(r(t))At/2 + Ar D ^At/2. The particle sits now at a 
position r(t + At) 

5. Consider the probability that a particle located at r(t + At) reacts during a time interval 
At/2. As in step 1, if the reaction takes place, set the reaction time to tR = t + At/2 



TABLE I: Algorithm for a coordinates dependent X — > process 
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Algorithm II: multiple time step algorithm 

1. Repeat J-/2 times 

l.a. If the particle is B, consider changing its state during a time interval At/2, else do nothing, 
l.b. If the particle is still B, proceed with algorithm I during the time interval At. 

1. e. If the particle is still B, consider changing its state during a time interval At/2, else do nothing. 

2. Do once 

2. a. If the particle is A, consider changing its state during a time interval TAt/2, else do nothing. 
2.b. If the particle is still A, proceed with algorithm I during the time interval TAt. 

I.e. If the particle is still A, consider changing its state during a time interval TAt/2, else do nothing. 

3. Repeat T /2 times 

l.a. If the particle is B, consider changing its state during a time interval At/2, else do nothing. 

l.b. If the particle is still B, proceed with algorithm I during the time interval At. 

I.e. If the particle is still B, consider changing its state during a time interval At/2, else do nothing. 



TABLE II: Multiple time step algorithm built on algorithm I 
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factor T 


time (s) 


zone [8 in f,6 sup ] 


A™ 


CPU time ratio 


1 


299.8 





— 


1.0 


4 


106.0 


[0.4-0.6] 


0.2 


0.354 


8 


62.0 


[0.4-0.6] 


0.2 


0.207 


12 


51.2 


[0.4-0.6] 


0.2 


0.171 


16 


42.6 


[0.4-0.6] 


0.2 


0.142 


16 


45.1 


[0.3-0.7] 


0.4 


0.150 


32 


33.5 


[0.3-0.7] 


0.4 


0.112 



TABLE III: Summary of CPU times relative to the absorption of 1000 independent particles at 6 C = 0.03. Brackets indicate 
the limits of the exchange zone. See text for details. 
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FIGURE CAPTIONS 
FIGURE 1. 

This picture shows the principle of a multiple time step procedure for a reaction-diffusion process. Finding the 
small reaction well requires the step size between two consecutive spatial positions to be smaller than the size of 
the well. Large steps might miss the reaction well completely, and the simulated reaction rate can only be largely 
underestimated. However, far from the reaction zone, large steps cause no harm, and can speed up the exploration 
of the configuration space by a large factor. Nevertheless, the boundary between the two zone (dashed circular line) 
will be subject to sampling artifacts if the procedure is done without care. 

FIGURE 2. 

This graph presents the spurious features expected from a naive multiple time step procedure in the absence of 
exchange zone. When an accelerating factor J- = 8 or J- = 32 is introduced, two peaks appear around x = ±0.2. 
Their width directly depend on the length step, namely \Z2DTAt. This artifact comes from a local breakdown of 
detailed balance induced by the time discretization. The special case T = 1, with no acceleration, gives a control 
curve, constant as expected, up to statistical fluctuations. 

FIGURE 3. 

Particle density in the absence of reaction when dynamics is governed by the internal state exchange dynamics in 
one dimension. "Slow" B particles occupy the central regions, and "fast" A particles mainly occupy the outside 
of the switching zone. The flat curves show the joined A and B distribution, which remains constant. When the 
exchange rates are increased, the A/ B proportion tends towards its limit value (jllj) . here a piecewise linear function. 
The exchange region where rates wab{t) and WBA{f) are both non zero is the interval [0.2,0.3] corresponding to an 
exchange zone size = 0.1. The acceleration factor used in this case is T = 8. The artifact in Fig. [2] is no longer 
visible. 
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FIGURE 4. 

Schematic view of the d s -dimensional configuration space, bounded by a reflecting outer boundary R and an absorbing 
inner boundary i?,. Two exchange zones are indicated by grey areas limited by dashed lines. 

FIGURE 5. 

Repartition (density) of particles as a function of r, with R e ~ 2.5, Ri = 1.0 and d s = 5. The exchange zone size A w 
here approaches 0.15. The inset shows an enlargement of the inner exchange zone. 

FIGURE 6. 

When the exchange zone is shrunk to A w = 0, under the same conditions as Fig. [5j the unphysical heterogeneity 
(artefact) appears near the outer exchange zone (r ~ 2.3) and, to a lesser extent, near the inner exchange zone 
(r ~ 1.15 as shown on the inset of this graph). 

FIGURE 7. 

Comparison of the survival rate, in Laplace space, between simulations and exact result, for spatial dimensions 
d s = 3,5,8 and 10 (see e.g. equations (|C6[) and (|C8|) in Appendix |C|) . We use for the simulations an acceleration 
factor T = 8 and a sample size of N = 1000 independent particles. The horizontal axis is the Laplace variable s while 
the vertical axis has the dimension of a time. The extrapolated value as s — > gives the average residence time (t) of 
the particles, and is subject to a statistical uncertainty of order (t) / \/~N (the distribution of residence times is nearly 
exponential) . 

FIGURE 8. 

Pictorial view of the reaction-diffusion problem on a sphere. The exchange zone A o B is indicated by the shaded 
area. Near the capture zone, B particles arc subject to the small step dynamics while apart from it, A particles are 
subject to the large step dynamics. 

FIGURE 9. 

The mean survival time (r r ) is calculated as a function of reaction zone size 9 C using different values of the accel- 
eration T (4 to 32 indicated as symbols). The simulation results show perfect agreements with analytic results: 
(r r ) = — ln(l — cos(0 c )) — 1 + In 2 + (1 — cos(# c ))/2, regardless of different values of J 7 . The time is measured in the 
unit of rotational diffusion coefficient D r Q. The exchange zone A w = 0.4 is chosen to be (6 C + 0.15, 8 C + 0.55). 
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FIGURE 10. 

Recurrence histogram of A and B particles as a function of angular positions 9. Particles disappear as soon as 
they encounter the capture radius 6 C = 0.03 and an exchange zone is included in the region 0.4 < 9 < 0.6. Solid 
lines correspond to A particle with T = 8, 16, 32 from top to bottom. Symbols are for B particles with the same 
acceleration factors. A stationary profile is expected to be reached when particles are injected as the same rate as 
they are captured. The inset displays possible artifact (kinks) near the exchange zone caused by too large T values, 
with fast particles A able to diffuse throughout the exchange zone and to reach 9 C . Nevertheless, the capture rate and 
the survival time distribution show little dependence on the choice of JF (cf Fig. [5]). 

FIGURE 11. 

Recurrence histogram of A and B particles in the interval of [9, 9 + d6] in the absence of reaction. The exchange zone 
is located in the interval 0.4 < 9 < 0.6. Bin counts are normalized by sin(0), so that uniform particle density on a 
sphere gives a flat histogram. The relative proportion of A and B particles changes smoothly across the exchange 
zone whereas the joined histogram remains flat. Examples are computed for T = 16 and T = 32. The total number 
of particles in the case of T — 16 is twice as large as that of T — 32, for given total exploration time. 

FIGURE 12. 

CPU time vs 1/F for values of F ranging between 4 and 256. It confirms the trend CPU time ~ a + (1 — ctj/J 7 with 
a fraction a ~ 0.06 (9 ~ 0.5), consistent with the exchange zone [0.4, 0.6]. Note that large values of T are not only 
time-inefficient, but also reintroduce the unwanted spurious features, calling for an optimal compromise. 

FIGURE 13. 

Reaction times are evaluated by three different numerical schemes: (*) a multiple-time-step procedure (J- = 16) 
without any exchange zone (naive approach), (o) a multiple-time-step procedure (F = 16) with exchange zone and 
single time step procedure T = 1 (dashed line). For the given number of receptors (Nr = 1, 5, 10), reaction times are 
measured function of the number of ligands Nl . 

FIGURE 14. 

Comparison of the computation times using three algorithms. In the presence of a single receptor Nr = 1, the 
computation times for 500 reactions are measured as a function of the number of ligands Nl ■ 

FIGURE 15. 

Comparison of computational times using three algorithms. For given Nr = 10, the computation time required for 
500 reactions are measured as a function of Nl. See text for details. 
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